function fig2_16

alpha=1;
lambda=1;
%b=sqrt(2)*12;
b=16;
t1=50;
t2=180;
nx=200;
x=linspace(-40,40,nx);
for ix=1:nx
    u0(ix)=exp(-x(ix)^2/b^2)*cos(lambda*x(ix));
end

x1=linspace(50,250,2*nx);
for ix=1:2*nx
    h1(ix)=sqrt(2/pi)*rhom2(x1(ix),t1,lambda,b);
end

x2=linspace(250,850,4*nx);
for ix=1:4*nx
    h2(ix)=sqrt(2/pi)*rhom2(x2(ix),t2,lambda,b);
end

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',0.01,'ML',0.052,'P',0.04)

hold on
plot(x,u0,'b','LineWidth',0.7)
plot(x1,h1,'r','LineWidth',0.7)
plot(x2,h2,'Color',[11/256, 218/256, 81/256],'LineWidth',0.8)
axis([-50 850 -1 1])
box on
grid on
%set(gca,'xtick',[0 1 2 3 4])

set(gca,'ytick',[-1 -0.5 0 0.5 1])
set(gca,'yTickLabel',{'$-1$';'$-0.5$';'0';'0.5';'1'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('$h(x,t)$','Interpreter','latex','FontSize',14)

legend('$\, t=0$','$\, t=50$','$\, t=180$','Location','NorthEast','FontSize',13,'Interpreter','latex','AutoUpdate','off')
plot([-50 900],[0 0],'k','LineWidth',0.5)
%exportgraphics(gcf,'/Users/mark/Desktop/kdv.pdf')

function RR=rhom2(x,t,lambda,b)
%  Integrate  over  [a, b]  using composite simpson and romberg
%  S(k) = composite simpson;  R(k) = romberg with composite simpson
%  err = relative iterative error using romberg

a=0; b=lambda+5/b;

%  interv = number of subinternals
tol=1e-8;
n=8;
S(1)=simp(a,b,n,x,t);
err=10*tol;
R(1)=1;  errR(1)=-1;
k=1;
while err>tol
    k=k+1;
    n=2*n;
    S(k)=simp(a,b,n,x,t);
    R(k)=(16*S(k)-S(k-1))/15;
    err=abs(1-R(k)/R(k-1));
end
RR=R(k);


function y=f(k,x,t)
lambda=1; alpha=1; b=16;
y=(((exp(-(k-lambda)^2*b^2/4)+exp(-(k+lambda)^2*b^2/4))*b*sqrt(2))/4)*cos(-alpha*k^3*t+k*x);


%%% composite simpson: n = number of subintervals
function s=simp(a,b,n,x,t)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=f(xd(1),x,t);
for j=2:2:n
    ff2=f(xd(j+1),x,t);
    s=s+4*f(xd(j),x,t)+2*ff2;
end
s=h*(s-ff2)/3;














